clear;
sizes = [0, 35, 50, 75, 100, 139, 150, 160, 200, 250, 300, 340, 350, 400, 450, 490, 500];

load result.txt;

% coorditanes of the ladder peaks on the original time axis
for i = 1:length(sizes)
  [d, x] = min(abs(X(1,:) - sizes(i)));
  xs(i) = x;
end
xs = [xs, length(X)]; % include the endpoint; 0 is already there (see sizes above)

% apply the warping path in result.txt to the sample and to the
% straight line, T
T = [1:length(X)];
for aN=1:diagnos.Nsegments-1
    tindex = diagnos.indexT(aN):diagnos.indexT(aN)+diagnos.segment_length-1;
    nt = length(tindex);
    pindex = warping(aN):warping(aN+1)-1;
    np = length(pindex);
    p = (0:(nt-1))*(np-1)/(nt-1) + 1;
    xw(tindex) = interp1(1:np,T(pindex),p);
end
tindex = diagnos.indexT(diagnos.Nsegments):length(X(2,:));
nt = length(tindex);
pindex = warping(diagnos.Nsegments):length(X(2,:));
np = length(pindex);
p = (0:(nt-1))*(np-1)/(nt-1) + 1;
xw(tindex) = interp1(1:np,T(pindex),p);

% calculate the warping differential (a difference between the
% straight line (1:length(xw)) and the warping path xw.
df = xw - T;

% The only interesting values of the differential are those
% at the ladder peaks. Interpolate everything in between
time = 1:length(X);
dfi = interp1(xs, df(xs), time);
dfi(end) = dfi(end-1); % get rid of the last Nan
warped_time = time - dfi;

% interpolate sample from warped to straight time
ws = interp1(warped_time, X(3,:), time);
ws(end) = ws(end-1);

% quadrature of the differential
for i = 1:length(df)
  q(i) = sum(df(1:i)); # quadrature
end

%plot it all
xseg = X(1,diagnos.indexT); % segment boundaries

subplot(2,1,1)
plot(
X(1,:), X(2,:), 'g;target;',
X(1,:), -X(3,:)-10, 'm;- raw sample;',
X(1,:), -Xw_cow(1:length(X)), 'b;- sapmle warped with COW;',
X(1,:), ws, 'r;sample warped with interpolated path;',
xseg, max(X(3,:))*0.05*ones(1,length(xseg)), '^k;segment boundaries;'
);


#axis([400,520]);
grid;
subplot(2,1,2)
plot(
X(1,:), df,'r;warping differential (straight line - warping path) [time points];',
X(1,:), dfi, 'm;interpolated warping differential [time points];',
X(1,:), 10*(df - dfi)*max(X(1,:))/xs(end-1), 'b; 10 x deviation [base pairs];',
X(1,xs), df(xs), '^*;target peak positions;',
X(1,:), 0*X(1,:), 'k;zero;'
);
#axis([400,520]);
ylabel("time points");
xlabel("length [base pairs]");
grid;

